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Abstract. We compute by numerical integration of the Dirac equation the 
number of quark-antiquark pairs produced in the classical color fields of collid- 
ing ultrarelativistic nuclei. Results for the dependence of the number of quarks 
on the strength of the background field, the quark mass and time are presented. 
We also perform several tests of our numerical method. While the number of 
qq pairs is parametrically suppressed in the coupling constant, we find that in 
this classical field model it could even be compatible with the thermal ratio to 
the number of gluons. 
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1. Introduction 



Due to large densities, implying large occupation numbers, the initial stages of an 
ultrarelativistic heavy ion collision may be be dominated by strong classical color 
fields. There is a twofold interest in calculating the production of quark-qntiquark 
pairs from these classical fields. Firstly, although heavy quark production is in 
principle calculable perturbatively, it would be interesting to understand whether 
these strong color fields influence the result. Secondly, being able to compute both 
gluon and quark production in the same framework would give insight into the 
chemical equilibration of the system and the consistency of the assumption of gluon 
dominance. The number of quark pairs present in the early stages of the system 
also has observable consequences in the thermal photon and dilepton spectrum. 
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Fig. 1. Domains of different time de- 
pendences. The fermion amplitude is a 
sum of two terms: one interacting first 
with the left moving nucleus, then the 
right moving one, and vice versa, ^i*' 2 ^ 
are pure gauges and A j , A v is a numeri- 
cally computed color field. In the pertur- 
bative limit of the Abelian theory these 
terms correspond to the u and t channel 
diagrams in Fig. 1121 

In this talk we shall present first results [ d 12] of a numerical computa- 
tion of quark antiquark pair production from the classical fields of the McLerran- 
Venugopalan (MV) model. The equivalent calculation, although in the covariant 
gauge unlike the present computation, has been carried out analytically to lowest 
order in the densities of both color sources ("pp"-case) in Ref. and to lowest 
order in one of the sources ("pA"-case) in Ref. The corresponding calculation 
in the Abelian theory [0 Ejj of interest for the physics of ultraperipheral collisions, 
can be done analytically to all orders in the electrical charge of the nuclei. Quark 
pair production has also been studied in a related "CGC" -approach in Refs. [d |5] 
and in a more general background field in Ref. 

2. The numerical calculation 

Our calculation of pair production relies on the numerical calculation of the classical 
background color field, in which we solve the Dirac equation. 

2.1. The background Geld 

In the classical field model the background gluon field is obtained from solving the 
Yang-Mills equation of motion with the classical color source J v given by transverse 
color charge distributions of the two nuclei boosted to infinite energy: 

[£V,FH = J" = S» + p {1) (x T )S(x-) + S' J -p ( 2)(xT)6(x + ). (1) 

In the MV model [^J the color charges are taken to be random variables with a 
Gaussian distribution 

(p a (x T )p b (y T )) = g 2 » 2 6 ab 6 2 (x T - y T ) (2) 

depending on the coupling g and a phenomenological parameter ji. The combination 
g 2 fi is closely related to the saturation scale Q s . Collisions of two ions were first 
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Fig. 2. Amplitude (in lattice units) at 
t = 0.25 fm as a function of Ay for differ- 
ent antiquark momenta on a 180 2 -lattice 
(see text in Sec. 12.41 for the notation). 



Fig. 3. Dependence on proper time r 
of the number of pairs of one flavor per 
unit rapidity dN/dy for g 2 \i = 2 GeV 
and different quark masses. The lowest 
curve corresponds to g 2 /j, — 1 GeV. 



studied analytically using this model in Ref. [ Ill| and the way of numerically solving 
the equations of motion in a Hamiltonian formalism in the A T = gauge was 
formulated in Ref. [IT2*). 



2.2. Solving the Dirac equation 

Our method of solving the Dirac equation is explained in more detail and the 
numerics tested in a 1+1-dimensional toy model in Ref. [ I13|. The domains of 
spacetime involved in the calculation are illustrated in Fig. ^ 

One starts in the the infinite past t — > — oo with a negative energy plane wave 
solution ip(x) = e lq ' x v(q). The Dirac equation can then be integrated forward in 
time analytically to the future light cone (r 2 = 2x + x~ = 0, x ± > 0) because the 
background field in the intermediate region is a pure gauge. This gives an initial 
condition, given explicitly in Eq. (16) of Ref. f° r ' l P{ T — 0, z,xr). Starting 

from this initial condition we then solve numerically the Dirac equation for r > 
using the coordinate system r, z, xy. The reason for this choice of coordinates is the 
following. It not feasible to have degrees of freedom at different energy scales (yfs 
and Q s ) present in the same numerical calculation. Thus the temporal coordinate 
is chosen to be r in order to include the hard sources of the color fields only in the 
initial condition. Although the background field is boost invariant, the longitudinal 
coordinate can not be disregarded in this computation because the rapidities of the 
quark and the antiquark are correlated. The longitudinal coordinate is chosen as z, 
not the usual dimensionless r\ = tanhT (z/t), because the initial condition on the 
light cone involves dimcnsionful longitudinal momentum scales (coming from the 
four momentum q), and in order to represent them on a spatial lattice at t = a 
dimensionful coordinate is needed. 

The pair production amplitude M T at proper time r is then obtained by fixing 
the Coulomb gauge in the transverse plane diAi — and projecting the spinor 
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Fig. 4. Dependence of the number of 
quark pairs on quark mass at a fixed 
proper time, r = 0.25 fm, and for two 
values of g 2 \i. 
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Fig. 5. Dependence of the number 
of quark pairs on g 2 \i at a fixed proper 
time, t = 0.25 fm, and for quark mass 
to = 0.3 GeV. 




wavefunction ip(r, z, x^) onto positive energy states 4> p {t, z, x^) = e lp ' x u(p): 

M T (p,q)= f ^ffiZ. 0t (t, z, x r )7V V«(t, z, x t ). (3) 

For times larger than the formation time of the quark pair r > 1/ \Jm 2 + q^ 2 this 
amplitude can be interpreted as the amplitude for producing quark antiquark pairs. 

2.3. Results 

Figure |21 shows \M T | 2 as a function of Ay = y p — y q integrated over pr for different 
q*r- When integrating also over the rapidity difference Ay one gets the number 
of pairs per unit rapidity as a function of r, shown in Fig. [3] It can be seen that 
the quark production amplitude reaches a finite value instantaneously and then 
increases slowly with r. 

The physical parameters of the calculation are g 2 /j characterizing the strength of 
the background field, the nuclear radius Ra and the quark mass to. The dependence 
on g 1 [i and to of the number of pairs at r = 0.25 fm is shown in Figs.21and|SJ The 
number of pairs is seen to increase with g 2 \i, but not as strongly as the (g 2 ^) 2 
dependence predicted by a simple dimensional analysis argument. The result also 
decreases with increasing quark mass, but the perturbative I/TO 2 -behavior is not 
reached in our calculation. The transverse momentum spectra of the (anti)quarks 
as a function of is shown for different quark masses and saturation scales in 
Figs. El and □ 

2.4. The numerical method 

Our numerical method is presented in more detail in Refs. The discretiza- 

tion in the transverse plane is straightforward, but in the longitudinal direction the 
r, z coordinate system can easily result in an unstable one. In Ref. [ I13j a stable 



Quark-antiquark production from classical fields 



5 




Fig. 6. Transverse momentum spec- 
trum of (anti)quarks for g 2 fi = 2 GeV 
at a fixed proper time, r = 0.25 fm, and 
for different quark masses. 



Fig. 7. Transverse momentum spec- 
trum of (anti)quarks for quark mass 
m = 0.3 GeV and for different g 2 \ju at 
a fixed proper time, r = 0.25 fm. 



scheme was found by using an implicit method, where at each timestep one solves 
a linear system of equations for the spinor at different points on the longitudinal 
lattice. 

The numerical computation depends on several discretization parameters: the 
lattice size and spacing in the longitudinal (N z and dz) and the transverse (N 2 and 
a) directions and the timestep dr. To check our numerical method we studied the 
dependence on these parameters. We have also tested the numerical method in the 
analytically solvable case of zero external field and tested how well our numerical 
implementation preserves boost invariance. 

The dependence on the longitudinal lattice parameters, dz and N z , is studied 
in Figs. |H1 and El Our convention is that the longitudinal lattice has points i = 
—N z . . . N z , i.e. 2N Z + 1 points and the length of the longitudinal lattice is 2N z dz. 
The results presented earlier in this talk are for dz = 0.2a and N z — 200 and have 
not been extrapolated to the infinite volume limit iV z dz — > oo. A finite dz leads to 
a lattice cutoff in p z and thus sets a maximum for the accessible interval in Ay. As 
can be seen from Fig. and understood as a consequence of the initial condition, 
Eq. (16) of Ref. [EH; the accessible region in Ay is smaller for small transverse 
momenta q-r- This cutoff in Ay is also apparent in Figs . 151 ITUl and ITTl 

We have so far only used transverse lattices of 180 2 points. The lattice momenta 
can be represented as = (q x ,q y ) with q x>y = —89. . .90, of which the g X)JI = 
—45 ... 45 are non-doubler modes. We explicitly leave out the doubler modes both 
in the initial condition (q^ modes) and in the projection to the positive energy 
state (jpt modes) . This limits the volume of the transverse momentum space to 1/4 
of the bosonic case. Also the spectrum of quarks (see Figs. and [7J decreases so 
slowly for large transverse momenta that some dependence on the 1/a lattice cutoff 
is expected. Further computations with larger transverse lattices are still needed to 
study this issue. 

The memory requirement for a 400 x 180 2 -lattice with 4 x (N c — 3) complex 
components in one spinor is 1.2 GB in single precision. This can still be achieved on 




Fig. 8. Amplitude at r = 0.25 fm for 
g 2 /d = 2 GeV for different values of dz 
and N z . The lattice spacing is mea- 
sured in units of a, the transverse lattice 
spacing. 



Fig. 9. Number of pairs as a function 
of (half) the volume of the longitudinal 
lattice N x dz for two different values of 
dz. The solid lines are fits of the form 
A + B/{N z dzf. 



one processor, but for larger lattices a parallelized version of the program will have 
to be used. The numerical computations have been done on the "ametisti" cluster, 
a 66 x 2 processor AMD Opteron Linux cluster at the University of Helsinki, using 
over 10 17 flop of CPU-time so far. 

As explained in Fig.^ the amplitude is a sum of two terms. For a zero external 
color field these two terms, with an absolute value \M T \ — 1/ cosh(Ay/2) have an 
opposite sign and cancel each other. Comparing the analytically and numerically 
computed amplitudes is a nontrivial test of the numerical computation. Figure 1101 
shows how the analytical result for one branch and the cancellation when both 
branches are included are reproduced by our numerical method for one value of qr 
(qr = (5, 5) on a 150 2 lattice). 

Due to the boost invariance of the background field the amplitude M T should 
not depend on the rapidities y q and y p separately, but only on the difference Ay = 
Vp — Dq- Because the calculation is done using z, not rapidity, as the longitudinal 
coordinate, verifying the boost invariance of the resulting amplitude is a nontrivial 
test of the numerical method. Figure ITT1 shows that the number of pairs is, taking 
into account numerical inaccuracies, not affected by the choice of y q . 



3. Discussion 



It has conventionally been assumed that the initial state of a heavy ion collision 
is dominated by gluons. This is the result e.g. when both quarks and gluons are 
produced in 2 — * 2 collisions of collinear partons [I15|. Indeed, when comparing the 
cross sections of the 2 — > 2 processes for quark pair production and gluon production 
(see Fig. El, the diagrams for quark pair production are suppressed by a factor of 
~ 210 — 7 x 30, of which the factor 7 is due to color factors. 

In the color glass condensate framework the picture is quite different. Whereas 
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Fig. 10. The amplitude in a zero exter- 
nal color field for two different timesteps 
dr. Plotted are the absolute value of 
the amplitude for only one and both 
branches of the amplitude (see text and 
Fig.©. 
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Fig. 11. The amplitude \M T \ 2 as a 
function of the rapidity difference Ay for 
different rapidities of the antiquark y q . 
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Fig. 12. The 2 — ► 2-processes for producing quark pairs and gluons from an initial 
gluon distribution. 

the quark pairs are produced in a 2 — * 2 process similarly to collinear factorisation, 
gluons are produced in what reduces in the perturbative limit to a 2 — > 1 process 
with a smaller power of the coupling g. It is therefore less straightforward to 
compare the two. Strictly perturbatively quark pair production is of higher order 
in g. But the phenomenologically relevant value is g ~ 2, and thus higher orders 
in g are not suppressed by much when looking at the actual numbers. One can 
then legitimately ask whether one should, for consistency, also include quantum 
corrections (higher orders in g) in the computation of gluon production. 

Brushing aside these remarks and boldly taking our numerical results at their 
face value would lead to the following scenario. It was earlier assumed that the initial 
state of a heavy ion collision is dominated by gluons. If the subsequent evolution of 
the system conserves entropy and thus, approximately, multiplicity this would mean 
~ 1000 gluons in a unit of rapidity. In the classical field model this corresponds [ 
116) to g 2 /i « 2 GeV. Our results seems to point to a rather large number of quark 
pairs present already in the initial state. One could envisage a scenario where for 
g 2 // w 1.3 GeV these 1000 particles could consist of > 400 gluons, > 300 quarks 
and > 300 antiquarks (take the lowest curve from Fig. [3] and multiply by Nf = 3). 
This would be close to the thermal ratio of N g /N q = 64/(21iV f ). 
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4. Conclusions 

We have calculated quark pair production from classical background field of McLerran- 
Venugopalan model by solving the 3+1-dimensional Dirac equation numerically in 
this classical background field. We find that number of quarks produced is large, 
pointing to a possible fast chemical equilibration of the system. The mass depen- 
dence of our result surprisingly weak and we are not yet able to make any conclusions 
on heavy quarks until studying the numerical issues involved. 
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